Fully reported case
Simulation
# q_constant_FR
params_q_constant_FR <- list(
data = list(
alpha_lamb = alpha_lamb,
beta_lamb = beta_lamb,
T = T,
date_start = first_date,
D = D
),
q_model = list(
method = "q_constant",
method_params = list(b = 0.7, phi = 0.9)
)
)
q_constant_FR <- simulateData(params_q_constant_FR)
plot(q_constant_FR$q, col = "red", pch = 19, type = "b", ylim = c(0,1))
Exploratory analysis
# exploritary analysis
page_num <- ceiling(nrow(q_constant_FR$case_reported_cumulated)/16)
exp_plot_q_constant <- fit_exp_plot(q_constant_FR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# print(exp_plot_q_constant)
#
# exp_plot_q_constant$coefficients
exp_b_data_q_constant<- data.frame( date = as.Date(rownames(q_constant_FR$case_reported_cumulated)),
b = exp_plot_q_constant$coefficients$b)
exp_b_plot_q_constant <- ggplot(exp_b_data_q_constant, aes(x = date, y = b)) +
geom_point(color = "black", size = 1.5) +
geom_smooth(method = "loess", se = TRUE,
color = "blue", fill = "grey", alpha = 0.5) +
theme_minimal() +
labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")
print(exp_b_plot_q_constant)
## `geom_smooth()` using formula = 'y ~ x'
Model fitting
# out_q_constant_FR <- nowcasting_moving_window(q_constant_FR$case_reported_cumulated, scoreRange = scoreRange,
# case_true = q_constant_FR$case_true,
# start_date = first_date,
# D = D,
# methods =models_to_use,
# compiled_models = compiled_models,
# iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
# num_chains = 3, suppress_output = T,
# posterior_draws_path = file.path(posterior_draws_path, "q_constant")
# )
#
#
# save(out_q_constant_FR, file = file.path(data_save_path, "FR_q_constant.RData"))
load( file.path(data_save_path, "FR_q_constant.RData"))
#
results_q_constant_FR <- nowcasts_table(out_q_constant_FR, D = D, report_unit = "day",
methods = models_to_use
)
plots_q_constant_FR <- nowcasts_plot(results_q_constant_FR, D = D, report_unit = "day",
methods = models_to_use
)
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
for (i in 1:length(results_q_constant_FR)) {
print(calculate_metrics(results_q_constant_FR[[i]],
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 93.28 66.31 59.47 64.84 146.71 0.6 q_constant
## 2 6.99 9.26 4.35 6.52 64.61 1.0 b_constant
## 3 7.53 9.32 4.57 6.53 80.10 1.0 b_rw
## 4 8.20 10.20 4.95 7.16 75.51 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 47.37 10.47 35.96 10.06 111.41 0.9 q_constant
## 2 18.34 2.95 8.61 2.15 87.40 1.0 b_constant
## 3 18.90 2.91 8.75 1.94 94.81 1.0 b_rw
## 4 19.40 2.93 8.38 1.90 93.16 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 22.75 2.74 16.15 2.54 113.54 1 q_constant
## 2 17.03 1.73 7.30 1.12 108.77 1 b_constant
## 3 19.65 1.81 7.73 1.20 118.37 1 b_rw
## 4 17.74 1.74 7.06 1.09 116.55 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 32.10 3.43 15.46 1.97 118.36 1 q_constant
## 2 25.47 2.78 7.96 1.17 115.94 1 b_constant
## 3 33.40 3.57 9.59 1.35 120.04 1 b_rw
## 4 25.10 2.77 7.71 1.13 117.81 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 15.03 3.09 9.59 1.66 110.35 1 q_constant
## 2 10.80 2.55 4.84 1.08 109.25 1 b_constant
## 3 11.57 2.79 4.64 1.08 109.91 1 b_rw
## 4 10.82 2.58 4.58 1.02 109.63 1 b_ou
for (i in 1:length(results_q_constant_FR)) {
print(calculate_metrics(data.table::last(results_q_constant_FR[[i]], D_check),
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 131.14 65.07 104.67 62.36 257.63 0.4 q_constant
## 2 9.85 12.59 8.01 9.93 110.01 1.0 b_constant
## 3 10.63 12.81 8.53 10.25 140.60 1.0 b_rw
## 4 11.57 13.98 9.20 11.05 131.61 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 71.70 10.22 67.86 9.53 244.84 0.8 q_constant
## 2 36.27 4.55 27.82 3.66 191.01 1.0 b_constant
## 3 37.38 4.68 28.57 3.74 220.83 1.0 b_rw
## 4 38.50 4.79 28.12 3.67 214.43 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 39.79 2.77 33.77 2.38 267.82 1 q_constant
## 2 40.61 2.66 30.06 2.07 257.42 1 b_constant
## 3 47.27 3.06 32.71 2.22 311.02 1 b_rw
## 4 42.74 2.77 30.75 2.09 295.44 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 84.52 8.62 56.44 5.65 231.21 1 q_constant
## 2 71.37 7.29 45.73 4.57 224.63 1 b_constant
## 3 93.93 9.61 58.47 5.87 256.82 1 b_rw
## 4 70.49 7.20 45.12 4.51 239.02 1 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 38.27 8.77 28.29 6.29 147.62 1 q_constant
## 2 32.00 7.46 22.47 5.12 147.22 1 b_constant
## 3 35.05 8.23 22.94 5.34 153.61 1 b_rw
## 4 32.60 7.60 22.65 5.17 149.21 1 b_ou
plots_q_constant_FR
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Parameter Checking
model q_constant
# try the third one, "2024-01-30"
T_test = T * 3/6
# q_constant
varnames_q_const <- out_q_constant_FR$q_constant[[3]]$summary()$variable
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_q_const, value = TRUE),
x = q_constant_FR$lambda[1:T_test]
)
mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^p\\[.+\\]$", varnames_q_const, value = TRUE),
x = c(q_constant_FR$q[1], diff(q_constant_FR$q))
)
mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("p"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[.+\\]$", varnames_q_const, value = TRUE),
x = q_constant_FR$q
)
mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("q"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_q_const, value = TRUE),
x = q_constant_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_FR$q_constant[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
model b_constant
# try the third one, "2024-01-30"
T_test = T * 3/6
# q_constant
varnames_b_const <- out_q_constant_FR$b_constant[[3]]$summary()$variable
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("b"), prob_outer = 0.95)
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("phi"), prob_outer = 0.95)
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_b_const, value = TRUE),
x = q_constant_FR$lambda[1:T_test]
)
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[.+\\]$", varnames_b_const, value = TRUE),
x = q_constant_FR$q
)
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("q"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_b_const, value = TRUE),
x = q_constant_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_FR$b_constant[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
Not fully reported
Simulation
# q_constant_NFR
params_q_constant_NFR <- list(
data = list(
alpha_lamb = alpha_lamb,
beta_lamb = beta_lamb,
T = T,
date_start = first_date,
D = D
),
q_model = list(
method = "q_constant",
method_params = list(b = 0.1, phi = 0.9)
)
)
q_constant_NFR <- simulateData(params_q_constant_NFR)
plot(q_constant_NFR$q, col = "red", pch = 19, type = "b", ylim = c(0,1))
Exploratory analysis
# exploritary analysis
page_num <- ceiling(nrow(q_constant_NFR$case_reported_cumulated)/16)
exp_plot_q_constant <- fit_exp_plot(q_constant_NFR$case_reported_cumulated,ncol = 4, nrow = 4, page = c(1:page_num), if_fit = T)
# print(exp_plot_q_constant)
#
# exp_plot_q_constant$coefficients
exp_b_data_q_constant<- data.frame( date = as.Date(rownames(q_constant_NFR$case_reported_cumulated)),
b = exp_plot_q_constant$coefficients$b)
exp_b_plot_q_constant <- ggplot(exp_b_data_q_constant, aes(x = date, y = b)) +
geom_point(color = "black", size = 1.5) +
geom_smooth(method = "loess", se = TRUE,
color = "blue", fill = "grey", alpha = 0.5) +
theme_minimal() +
labs(x = NULL, y = "Y", title = "Smoothed Curve of parameter b")
print(exp_b_plot_q_constant)
## `geom_smooth()` using formula = 'y ~ x'
Model fitting
# out_q_constant_NFR <- nowcasting_moving_window(q_constant_NFR$case_reported_cumulated, scoreRange = scoreRange,
# case_true = q_constant_NFR$case_true,
# start_date = first_date,
# D = D,
# methods =models_to_use,
# compiled_models = compiled_models,
# iter_sampling = 2000, iter_warmup = 1000, refresh = 0,
# num_chains = 3, suppress_output = T,
# posterior_draws_path = file.path(posterior_draws_path, "q_constant")
# )
#
#
# save(out_q_constant_NFR, file = file.path(data_save_path, "NFR_q_constant.RData"))
load( file.path(data_save_path, "NFR_q_constant.RData"))
results_q_constant_NFR <- nowcasts_table(out_q_constant_NFR, D = D, report_unit = "day",
methods = models_to_use
)
plots_q_constant_NFR <- nowcasts_plot(results_q_constant_NFR, D = D, report_unit = "day",
methods = models_to_use
)
for (i in 1:length(results_q_constant_NFR)) {
print(calculate_metrics(results_q_constant_NFR[[i]],
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 24.31 24.80 18.22 22.68 94.21 1.0 q_constant
## 2 66.30 41.02 44.27 37.72 57.50 0.5 b_constant
## 3 72.26 43.68 48.16 40.15 61.00 0.5 b_rw
## 4 70.92 43.60 47.44 40.36 59.61 0.5 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 120.02 24.24 89.48 22.20 98.36 0.30 q_constant
## 2 76.18 16.46 54.25 15.28 118.32 0.75 b_constant
## 3 96.02 18.39 66.22 16.76 157.52 0.80 b_rw
## 4 82.14 17.00 58.28 15.83 146.62 0.85 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 187.91 21.64 144.55 20.47 108.40 0.20 q_constant
## 2 38.54 11.16 25.60 7.25 138.21 0.97 b_constant
## 3 44.87 9.44 29.19 7.00 243.13 1.00 b_rw
## 4 37.40 10.53 24.44 6.79 193.75 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 193.25 20.41 161.35 19.56 111.08 0.17 q_constant
## 2 37.82 10.28 26.93 6.32 137.36 1.00 b_constant
## 3 30.79 8.30 22.80 5.52 212.57 1.00 b_rw
## 4 34.59 9.96 23.56 5.86 178.71 1.00 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 185.44 20.75 157.32 19.94 102.59 0.16 q_constant
## 2 27.47 9.34 20.93 5.65 121.27 1.00 b_constant
## 3 34.51 8.85 24.17 5.85 155.83 0.98 b_rw
## 4 25.23 8.90 19.15 5.26 149.45 1.00 b_ou
for (i in 1:length(results_q_constant_NFR)) {
print(calculate_metrics(data.table::last(results_q_constant_NFR[[i]],D_check),
methods = models_to_use))
}
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 33.89 18.54 31.16 18.31 162.42 1.0 q_constant
## 2 93.52 48.32 83.17 47.53 97.21 0.2 b_constant
## 3 101.96 52.51 90.85 51.84 104.01 0.2 b_rw
## 4 100.05 51.67 89.15 50.97 101.01 0.2 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 199.32 25.71 193.01 25.06 213.02 0.0 q_constant
## 2 128.61 16.60 117.95 15.28 258.24 0.8 b_constant
## 3 170.92 21.76 158.90 20.43 391.83 0.8 b_rw
## 4 140.98 18.11 130.80 16.90 343.64 0.8 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 350.03 24.54 344.44 24.28 251.01 0.0 q_constant
## 2 83.81 6.02 76.26 5.42 331.42 0.8 b_constant
## 3 93.44 6.67 86.99 6.17 705.27 1.0 b_rw
## 4 81.70 5.82 75.61 5.35 502.64 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 198.41 18.45 192.79 17.91 215.42 0.2 q_constant
## 2 79.50 7.73 71.94 6.84 276.82 1.0 b_constant
## 3 55.72 5.59 50.68 4.93 513.46 1.0 b_rw
## 4 71.23 6.72 59.57 5.61 393.62 1.0 b_ou
## RMSE RMSPE MAE MAPE Interval_Width Coverage_Rate Method
## 1 114.11 22.18 106.33 20.91 140.42 0.4 q_constant
## 2 43.82 9.58 37.29 7.92 170.21 1.0 b_constant
## 3 59.74 12.24 55.14 11.16 233.23 1.0 b_rw
## 4 41.26 8.94 35.22 7.42 201.01 1.0 b_ou
plots_q_constant_NFR
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
Parameter Checking
model q_constant
# try the third one, "2024-01-30"
T_test = T * 3/6
# q_constant
varnames_q_const <- out_q_constant_NFR$q_constant[[3]]$summary()$variable
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_q_const, value = TRUE),
x = q_constant_NFR$lambda[1:T_test]
)
mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
param_true = tibble(
parameter = grep("^p\\[.+\\]$", varnames_q_const, value = TRUE),
x = c(q_constant_NFR$q[1], diff(q_constant_NFR$q))
)
mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("p"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[.+\\]$", varnames_q_const, value = TRUE),
x = q_constant_NFR$q
)
mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("q"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_q_const, value = TRUE),
x = q_constant_NFR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_NFR$q_constant[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
model b_constant
# try the third one, "2024-01-30"
T_test = T * 3/6
# q_constant
varnames_b_const <- out_q_constant_NFR$b_constant[[3]]$summary()$variable
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("b"), prob_outer = 0.95)
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("phi"), prob_outer = 0.95)
param_true = tibble(
parameter = grep("^lambda\\[.+\\]$", varnames_b_const, value = TRUE),
x = q_constant_NFR$lambda[1:T_test]
)
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("lambda"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^q\\[.+\\]$", varnames_b_const, value = TRUE),
x = q_constant_NFR$q
)
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("q"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
param_true = tibble(
parameter = grep("^N\\[.+\\]$", varnames_b_const, value = TRUE),
x = q_constant_FR$case_true[1:T_test, 1]
)
mcmc_areas(out_q_constant_NFR$b_constant[[3]]$draws("N"), prob_outer = 0.95) +
geom_point(aes(x = x), param_true, color = "red", size = 1)
